pred_mnl <- function(formula, 
                      data, 
                      weights = NULL,
                      imp = FALSE,
                      S = 1000L,
                      seed = 19890213,
                      reorder.y = FALSE,
                      x.name,
                      x.steps = NULL,
                      x.shift = "unit",
                      ey = TRUE,
                      ey_stack = TRUE,
                      dydx = TRUE,
                      rdydx = FALSE) { 
  ## Seed
  set.seed(seed)
  
  ## Model Estimation
  print("1 Estimating models")
  if (imp) {
    M <- length(data)
    b <- v <- list()
    for (i in seq_len(M)) {
      df <- data[[i]]
      w <- if (is.null(weights)) {
        NULL
      } else {
        df[, weights]
      }
      df$weights <- w
      m <- multinom(formula,
                    Hess = TRUE,
                    data = df,
                    weights = weights)
      v[[i]] <- vcov(m)
      b[[i]] <- coef(m)
    }
    rm(df)
  } else {
    M <- 1L
    w <- if (is.null(weights)) {
      NULL
    } else {
      data[, weights]
    }
    data$weights <- w
    m <- multinom(formula,
                  Hess = TRUE,
                  data = data,
                  weights = weights)
    b <- coef(m)
    v <- vcov(m)
  }
  
  ## Extract model information
  call <- as.formula(paste(as.character(m$terms)[c(1, 3)], sep = " "))
  j.names <- m$lab
  k.names <- m$coefnames
  weights <- m$weights
  
  ## Constants
  J <- nrow(coef(m)) + 1L
  K <- ncol(coef(m))
  
  ## Parameter Simulation
  print("2 Simulating Model Parameters")
  if (imp) {
    b.sim <- NULL
    for(i in seq_len(M)) {
      coef.vec <- as.vector(t(b[[i]]))
      b.sim.tmp <- mvrnorm(S, coef.vec, v[[i]])
      for (k in seq_len(K)) b.sim.tmp <- cbind(0, b.sim.tmp)
      b.sim <- rbind(b.sim, b.sim.tmp)
      rm(b.sim.tmp)
      b.est <- rbind(0, apply(simplify2array(b), 1:2, mean))
    }
  } else {
    coef.vec <- as.vector(t(b))
    b.sim <- mvrnorm(S, coef.vec, v)
    for (k in seq_len(K)) b.sim <- cbind(0, b.sim)
    b.est <- rbind(0, b)
  }
  
  ## Simulations
  b <- array(NA, c(J, M * S, K))
  for (j in seq_len(J)) {
    from <- (j - 1L) * K + 1L
    to <- j * K
    b[j, , ] <- b.sim[, from:to]
  }
  rm(b.sim)
  
  ## Data matrices
  if (imp) {
    X <- lapply(data, function(d) model.matrix(call, data = d))
    if (is.null(x.name)) {
      X <- X[[M]]
    } else {
      X <- simplify2array(X) 
    }
  } else {
    X <- model.matrix(call, data = data)
  }
  N <- nrow(X)
  
  ## Find x-variables of interest
  if (is.null(x.name)) {
    n.cat <- 0L
    x.type = "bin"
  } else {
    if (M > 1) {
      X.tmp <- X[, , 1]
    } else {
      X.tmp <- X
    }
    x.index <- grepl(x.name, colnames(X.tmp))
    x.which <- which(x.index)
    if (sum(x.index) > 1) {
      x.type <- "cat"
      n.cat <- sum(x.index)
    } else if (length(unique(X.tmp[, x.index])) > 2) {
      x.type <- "cont"
    } else {
      x.type <- "bin"
      n.cat <- 1L
    }
    rm(X.tmp)
  }
  
  ## Average X over imputations
  if (M > 1 & !is.null(x.name)) X <- apply(X, 1:2, mean)
  
  ## Predictions
  print("3 Calculating requested quantities")
  if (x.type %in% c("bin", "cat")) { # discrete/no predictor
    ## Data
    X.tmp <- X
    if (!is.null(x.name)) {
      X.fill <- matrix(0, N, n.cat)
      X.tmp[, x.which] <- X.fill
      X.tmp <- replicate(n.cat + 1, X.tmp)
      for (i in 1:n.cat) {
        X.tmp[, x.which[i], i + 1] <- 1
      }
    } else {
      X.tmp <- replicate(n.cat + 1, X.tmp)
    }
    
    ## Probabilities
    if (ey | ey_stack | dydx) {
      pr.sim <- exp.xb.sim <-  array(NA, c(J, M * S, N, n.cat + 1))
      pr.est <- exp.xb.est <-  array(NA, c(J, N, n.cat + 1))
      for (c in seq_len(n.cat + 1)) {
        for (j in seq_len(J)) {
          exp.xb.sim[j, , , c] <- exp(b[j, , ] %*% t(X.tmp[, , c]))
          exp.xb.est[j, , c] <- exp(X.tmp[, , c] %*% as.matrix(b.est[j, ]))
        }
        for (j in seq_len(J)) {
          denom <- apply(exp.xb.sim[, , , c], 2:3, sum)
          pr.sim[j, , , c] <- exp.xb.sim[j, , , c] / denom
          denom <- apply(exp.xb.est[, , c], 2, sum)
          pr.est[j, , c] <- exp.xb.est[j, , c] / denom
        }
      }
      pr.sim <- apply(pr.sim, c(1, 2, 4), weighted.mean, w)
      pr.est <- apply(pr.est, c(1, 3), weighted.mean, w)
      
      pr.yx <- array(NA, c(3L, J, n.cat + 1))
      pr.yx[1, , ]   <- pr.est
      pr.yx[2:3, ,]  <-  apply(pr.sim, c(1, 3), quantile, c(.025, .975)) 
      
      ## Labeling
      if(is.null(x.name)) {
        dimnames(pr.yx) <-  
          list(c("Est.", "2.5%", "97.5%"),
               j.names)
      } else {
        dimnames(pr.yx) <-
          list(c("Est.", "2.5%", "97.5%"),
               j.names, 
               paste(x.name, seq_len(n.cat + 1), sep = " = ")) 
      }
    } else {
      pr.yx <- NULL
    }  
    
    ## Stacked Probabilities
    if (ey_stack) {
      pr.stk <- array(NA, dim(pr.yx))
      pr.stk[, 1, ] <- 0
      tmp <- 0
      for (j in 2:J) {
        tmp <- pr.sim[j - 1, , ] + tmp
        pr.stk[1, j, ]   <- pr.yx[1, j - 1, ] + pr.stk[1, j - 1, ]
        if (is.null(x.name)) {
          pr.stk[2:3, j, ] <- quantile(tmp, c(.025, .975))     
        } else {
          pr.stk[2:3, j, ] <- apply(tmp, 2, quantile, c(.025, .975))          
        }
      }
      rm(tmp)
      
      ## Labeling
      if(is.null(x.name)) {
        dimnames(pr.stk) <-  
          list(c("Est.", "2.5%", "97.5%"),
               j.names)
      } else {
        dimnames(pr.stk) <-
          list(c("Est.", "2.5%", "97.5%"),
               j.names, 
               paste(x.name, seq_len(n.cat + 1), sep = " = ")) 
      }
    } else {
      pr.stk <- NULL
    }
    
    ## Marginal Effects
    if (dydx & !is.null(x.name)) {
      dy.dx <- array(NA, dim = c(3L, J, n.cat))
      for (c in seq_len(n.cat)) {
        dy.dx[1, , c]   <- pr.yx[1, , c + 1] - pr.yx[1, , c]
        dy.dx[2:3, , c] <- apply(pr.sim[, , c + 1] - pr.sim[, , 1], 
                                 1, quantile, c(.025, .975))
      } 
      
      ## Labeling
      dimnames(dy.dx) <- 
        list(c("Est.", "2.5%", "97.5%"),
             j.names, 
             paste(x.name, 2:(n.cat + 1), sep = ": 1 vs "))
    } else {
      dy.dx <- NULL
    }
    
  } else { # continuous predictor
    
    ## Probabilities
    if (ey | ey_stack) {
      pr.sim <- exp.xb.sim <-  array(NA, c(J, M * S, N, length(x.steps)))
      pr.est <- exp.xb.est <-  array(NA, c(J, N, length(x.steps)))
      for (c in seq_along(x.steps)) {
        X.tmp <- X
        X.tmp[, x.which] <- x.steps[c]
        for (j in seq_len(J)) {
          exp.xb.sim[j, , , c] <- exp(b[j, , ] %*% t(X.tmp))
          exp.xb.est[j, , c] <- exp(X.tmp %*% as.matrix(b.est[j, ]))
        }
        for (j in seq_len(J)) {
          denom <- apply(exp.xb.sim[, , , c], 2:3, sum)
          pr.sim[j, , , c] <- exp.xb.sim[j, , , c] / denom
          denom <- apply(exp.xb.est[, , c], 2, sum)
          pr.est[j, , c] <- exp.xb.est[j, , c] / denom
        }
      }
      pr.sim <- apply(pr.sim, c(1, 2, 4), weighted.mean, w)
      pr.est <- apply(pr.est, c(1, 3), weighted.mean, w)
      
      pr.yx <- array(NA, c(3L, J, length(x.steps)))
      pr.yx[1, , ]   <- pr.est
      pr.yx[2:3, ,]  <-  apply(pr.sim, c(1, 3), quantile, c(.025, .975)) 
      
      ## Labeling
      dimnames(pr.yx) <- 
        list(c("Est.", "2.5%", "97.5%"),
             j.names, 
             paste(x.name, x.steps, sep = " = "))
    } else {
      pr.yx <- NULL
    }
    
    ## Stacked Probabilities
    if (ey_stack) {
      pr.stk <- array(NA, dim(pr.yx))
      pr.stk[, 1, ] <- 0
      tmp <- 0
      for (j in 2:J) {
        tmp <- pr.sim[j - 1, , ] + tmp
        pr.stk[1, j, ]   <- pr.yx[1, j - 1, ] + pr.stk[1, j - 1, ]
        pr.stk[2:3, j, ] <- apply(tmp, 2, quantile, c(.025, .975))
      }
      rm(tmp)
      
      ## Labeling
      dimnames(pr.stk) <- 
        list(c("Est.", "2.5%", "97.5%"),
             j.names, 
             paste(x.name, x.steps, sep = " = "))
    } else {
      pr.stk <- NULL
    }
    
    ## Marginal Effects
    if (dydx | rdydx) {
      ## Probabilities
      pr.sim <- exp.xb.sim <-  array(NA, c(J, M * S, N, 2L))
      pr.est <- exp.xb.est <-  array(NA, c(J, N, 2L))
      for (c in 1:2) {
        X.tmp <- X
        if (x.shift == "unit") {
          shift <- 1L
        } else if (x.shift == "sd") {
          shift <- sd(X.tmp[, x.which])
        } else if (x.shift == "tiny") {
          shift <- sd(X.tmp[, x.which]) / 1000
        }
        X.tmp[, x.which] <- X.tmp[, x.which] + (c - 1.5) * shift
        for (j in seq_len(J)) {
          exp.xb.sim[j, , , c] <- exp(b[j, , ] %*% t(X.tmp))
          exp.xb.est[j, , c] <- exp(X.tmp %*% as.matrix(b.est[j, ]))
        }
        for (j in seq_len(J)) {
          denom <- apply(exp.xb.sim[, , , c], 2:3, sum)
          pr.sim[j, , , c] <- exp.xb.sim[j, , , c] / denom
          denom <- apply(exp.xb.est[, , c], 2, sum)
          pr.est[j, , c] <- exp.xb.est[j, , c] / denom
        }
      }
      pr.sim <- apply(pr.sim, c(1, 2, 4), weighted.mean, w)
      pr.est <- apply(pr.est, c(1, 3), weighted.mean, w)
      
      dy.dx <- array(NA, dim = c(3L, J))
      dy.dx[1, ]   <- (pr.est[, 2] - pr.est[, 1]) / abs(shift)
      dy.dx[2:3, ] <- apply((pr.sim[, , 2] - pr.sim[, , 1]) / abs(shift), 
                            1, quantile, c(.025, .975))
      
      ## Labeling
      dimnames(dy.dx) <- 
        list(c("Est.", "2.5%", "97.5%"),
             j.names)
    } else {
      dy.dx <- NULL
    }
    
    ## Relative Marginal Effects
    if (rdydx) {
      div.sim <- apply(pr.sim, c(1, 2), mean)
      div.est <- apply(pr.est, 1, mean)
      
      rdy.dx <- array(NA, dim = c(3L, J))
      rdy.dx[1, ]   <- 
        ((pr.est[, 2] - pr.est[, 1]) / div.est) / (abs(shift))
      rdy.dx[2:3, ] <- 
        apply(((pr.sim[, , 2] - pr.sim[, , 1]) / div.sim) /abs(shift),
              1, quantile, c(.025, .975))
      
      ## Labeling
      dimnames(rdy.dx) <- 
        list(c("Est.", "2.5%", "97.5%"),
             j.names)
    } else {
      rdy.dx <- NULL
    }
  }
  
  ## Value
  print("4 Returning Output")
  out <- list(pr.yx  = pr.yx,
              pr.stk = pr.stk,
              dy.dx  = dy.dx,
              rdy.dx = rdy.dx)
  return(out)
}